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Characterization of the Earth’s surface is crucial to remote sensing, both to map 
geomorphological features and because subtracting this signal is essential during 
retrievals of the atmospheric constituents located between the surface and the sensor. 
Current operational algorithms model the surface total reflectance through a weighted 
linear combination of a few geometry-dependent kernels, each devised to describe a 
particular scattering mechanism. The information content of these measurements is 
overwhelmed by that of instruments with polarization capabilities: proposed models in 
this case are based on the Fresnel reflectance of an isotropic distribution of facets. 
Because of its remarkable lack of spectral contrast, the polarized reflectance of land 
surfaces in the shortwave infrared spectral region, where atmospheric scattering is 
minimal, can be used to model the surface also at shorter wavelengths, where aerosol 
retrievals are attempted based on well-established scattering theories. 

In radiative transfer simulations, straightforward separation of the surface and 
atmospheric contributions is not possible without approximations because of the 
coupling introduced by multiple reflections. Within a general inversion framework, 
the problem can be eliminated by linearizing the radiative transfer calculation, and 
making the Jacobian (i.e., the derivative expressing the sensitivity of the reflectance 
with respect to model parameters) available at output. We present a general metho- 
dology based on a Gauss-Newton iterative search, which automates this procedure and 
eliminates de facto the need of an ad hoc atmospheric correction. 

In this case study we analyze the color variations in the polarized reflectance 
measured by the NASA Goddard Institute of Space Studies Research Scanning Polari- 
meter during a survey of late-season snowfields in the High Sierra. This insofar unique 
dataset presents challenges linked to the rugged topography associated with the alpine 
environment and a likely high water content due to melting. The analysis benefits from 
ancillary information provided by the NASA Langley High Spectral Resolution Lidar 
deployed on the same aircraft. 

The results obtained from the iterative scheme are contrasted against the surface 
polarized reflectance obtained ignoring multiple reflections, via the simplistic subtrac- 
tion of the atmospheric scattering contribution. Finally, the retrieved reflectance is 
modeled after the scattering properties of a dense collection of ice crystals at the 
surface. Confirming that the polarized reflectance of snow is spectrally flat would allow 
to extend the techniques already in use for polarimetric retrievals of aerosol properties 
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over land to the large portion of snow-covered pixels plaguing orbital and suborbital 
observations. 

© 2012 Elsevier Ltd. All rights reserved. 


1. Introduction 

Meaningful retrievals of surface and atmospheric 
properties from space- or air-borne remote sensing data 
rely on accurate radiative transfer models as well as on 
successful separation of the surface and the atmospheric 
signals. In the case of retrievals over land, most opera- 
tional algorithms currently in use model the surface 
contribution by decomposing the total surface reflectance 
R in a sum of weighted kernels 

K 

R (8 s ,B r ,(j),b) = 'EhWMrJ, b). (1) 

k = 1 

The kernels I( k depend on the geometry of observation 
(i.e., on the solar zenith angle 9 S , the viewing zenith angle 
d r , and the Sun-sensor relative azimuth cp). Following 
the notation of Rodgers [1], the term b includes all other 
parameters involved non-linearly, such as wavelength 
(via the refractive index of the underlying medium) or 
the aspect ratio of the scattering objects at the surface. 
These dependencies, as well as that of the kernel weights 
f k on the wavelength, are dropped from the notation in 
what follows but remain understood. Although several 
kernels have been proposed, a common combination is 
the RossThick-LiSparse form 

k — f i s0 “F f vo l Kv„l “F f geo^geo (2) 

The first term on the right side of the equation is a 
Lambertian contribution with no directional dependence 
(i.e., kernel-free). The second term involves a kernel 
representing scattering within dense vegetation canopies 
and was derived by Roujean [2] from previous work [3], 
The last, geometric kernel accounts for shadowing of a 
Lambertian background from objects of larger size, and 
dates back to the model of Li and Strahler [4], although 
slightly modified by Lucht et al. [5], The functional form of 
these kernels can be found in the literature [6,7], 

The RossThick-LiSparse combination is used in retrie- 
val algorithms implemented for Earth-orbiting remote 
sensors such as the MODerate resolution Imaging Spec- 
troradiometer (MODIS [8,9,5]) and the POLarization and 
Directionality of the Earth’s Reflectances (POLDER [10]). 
The Multi-angle Imaging SpectroRadiometer (MISR) uses 
a different parameterization named Modified Rahman- 
Pinty-Verstraete (MRPV) model [11]. Additional Fresnel- 
type kernels [12,13] have been developed for POLDER, to 
this date the only spaceborne platform measuring polar- 
ization in the visible and near-infrared portion of the 
electromagnetic spectrum. The same kernels are also 
an option for modeling the data collected by the NASA 
Goddard Institute for Space Studies (GISS) Research Scan- 
ning Polarimeter (RSP) [14,15], the airborne prototype of 
the Aerosol Polarimetry Sensor (APS) which unfortunately 


failed to reach orbit due to a rocket underperformance at 
launch in 2011. 

The capability of measuring polarization adds greatly 
to the information content of remote sensing observa- 
tions, and can help to resolve ambiguities in retrievals 
that use intensity only both over land and over ocean 
[15-19]. Exploiting this type of information requires the 
understanding of the physical processes which can alter 
the essentially unpolarized state of the solar light incident 
at the top of the atmosphere. 

The information on the linear state of polarization of an 
observed scene is contained in the Q and U components of 
the Stokes vector. These two unbounded quantities are 
always specified relative to a reference system. Single 
scattering properties are naturally described relative to 
the scattering plane, which contains the directions of the 
incident and the scattered beam. As a simple consequence 
of the form of the phase matrix [20], the U component in 
this situation is inherently null so that the Q_ parameter 
alone is sufficient to specify the polarized reflectance 

R p = ^Jq^ + u 2 . For an instrument like the RSP, the implicit 
reference system is instead the meridional (scanning) 
plane. Even in this case, the Q parameter can be sufficient 
provided a rotation from the meridional into the scattering 
plane is first applied [21,16], However, for airborne obser- 
vations, an accurate estimate of the angle of rotation can 
suffer from uncertainties associated with instabilities in 
aircraft attitude. More importantly, the total observed 
signal is the sum of direct-beam and diffuse contributions: 
in problems involving surface reflection and multiple 
scattering in general a unique scattering plane for both 
these components cannot be defined unless the observa- 
tions take place along the principal plane, defined as the 
plane containing the incident direct beam and the local 
normal to the surface. At any other viewing geometry, 
rotation into either reference system leaves a residual U 
component. This contribution becomes detectable at view- 
ing angles close to backscatter, where both molecular and 
surface contributions are small and important structural 
signatures determined by aerosol (multiple) scattering 
may occur. The polarized reflectance has the modeling 
advantage of being independent of the reference system 
and of preserving the information contained in the resi- 
dual U component. It is nevertheless a positive quantity 
which often exhibits discontinuous derivatives at zero (in 
correspondence of polarization inversions), where care 
must be used in the statistical interpretation of the 
error bars. 

For most viewing geometries, the largest contribution 
to the polarized reflectance observed in the visible part of 
the electromagnetic spectrum from a downward-looking 
remote sensor comes from molecules and aerosols in the 
atmosphere. Rayleigh scattering from molecules produces 
light which vibrates perpendicularly to the plane of 
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scattering, while aerosols can polarize light both in the 
parallel and the perpendicular direction. The contribution 
due to the Earth's surface is smaller because it originates 
from reflection of the diffuse downward component, 
unless the observation takes place close to the principal 
plane. In a cone centered around the specular reflection 
direction, the surface can in fact have an overwhelming 
effect because of the high polarization introduced by the 
reflection of the direct beam. This condition is evident for 
sunglint over the ocean [22-24], governed by Fresnel 
reflection from a collection of facets typically described 
by a wave slope distribution [25], 

Specular reflection is also the basic mechanism gov- 
erning the reflection properties of most land surfaces, a 
fact that has been attributed to the smooth wax coating 
the leaves [26,27] and to the mineral facets of bare soils 
[28]. In a reflection event, the Fresnel formulae predict 
that the component of light perpendicular to the plane of 
reflection always exhibits higher reflectivity than does the 
parallel component, in a variable amount that depends on 
the angle of incidence. For this reason, a specular reflector 
always polarizes light in a direction perpendicular to the 
plane of reflection, with the extreme case of 100% of 
polarization introduced at Brewster angle. The fractional 
polarized reflectance F p (0,) is [29] 


F P (Vi) 


1 

2 


( rq cos 0i~n 2 cos 9 t \ 
\n ! cos 0j + n 2 cos 9 t J 


f n 2 cos 0,-rq cos 9 t \ 2 
\n 2 cos 0, + rq cos 9 t ) 


(3) 


shorter wavelengths. This approach is at the heart of 
techniques for the retrieval of aerosol properties over land 
based on high-accuracy polarimetric measurements [15], 

In summary, an appropriate polarized reflectance 
model should be able to closely mimic the surface signal 
with a twofold purpose: (a) to infer information on the 
surface characteristics such as vegetation cover [32], leaf 
orientation [33] and size [34], or macroscopic roughness 
[35]; and (b) to correct the signal when attempting to 
retrieve properties of the atmospheric constituents above 
the surface itself [36,15], 

This investigation tests the performance of the 
RossThick-LiSparse model applied to the unique surface 
type represented by alpine snow. The main goal, however, 
is to provide information on the spectral variability of the 
polarized reflectance of the surveyed snow surfaces over 
the 410-2264 nm range. Section 2 presents the iterative 
scheme used to retrieve an estimate of the surface signal, 
and its application to a synthetic scene. The following 
sections contain a description of the research flight 
performed over a snow target in June 2010 (Section 3), 
and a general summary of the optical properties of snow 
(Section 4). The results obtained by applying the iterative 
scheme to the collected dataset are contrasted with a 
standard atmospheric correction in Section 5. The conclu- 
sions are found in Section 6. 

2. Iterative atmospheric correction scheme 


where the angle of refraction 6 t is connected to the angle 
of incidence 0> by the Snell’s law of reflection 

jq sin = n 2 cos d t (4) 

and i q and n 2 are the indices of refraction of the 
two media. 

Corrections to the simple Fresnel model over land 
account for geometric effects such as multiple reflections 
and self-shadowing, so that ultimately some function of 
the Fresnel reflectance is used. Proposed models are based 
on the assumption of an isotropic distribution of facets on 
surfaces dominated by vegetation or bare soil [30,28,26]: 

R ,,e s= Fp(gi) • (5) 

p 4[cos(0 s ) + cos(0„)] ’ 


Fpjdj) 

4 cos(9 s )cos(6 v )' 


( 6 ) 


The independent variables of the model function (2) 
are the M values of the kernels I<j for the particular 
geometry of observation i. If the surface reflectance R s i 
were directly measured by the sensor, with hypothetical 
measurement error e f , it could be modeled as 

F s ,i = P(f + .Rj,i) + ti (7) 

or in matrix form, for JV viewing geometries as is the case 
for a scanning instrument: 

R s = K T f t + e (8) 

where f + is the Mxl state vector composed of an optimal 
choice of the weighting parameters f for each K,, now 
column vectors each of dimensions JV grouped into the 
I VI xN matrix K. Since the model is linear in the variables 
Kj, the solution could be found with a linear least-square 
regression formally expressed through normal equations 


The functional form for the polarized reflectance of the 
soil model of Eq. (6) has the disadvantage of leading to 
unphysical values of albedo (greater than 1) at some 
viewing geometries. In order to fit the surface polarized 
reflectance we employ the vegetation model, which 
instead ensures proper behavior, scaled with the weight 
f fresnei to account for variability in refractive index to which 
Eq. (5) is roughly proportional. 

As stated above, the real part of the index of refraction of 
most natural surfaces varies very little across the visible and 
shortwave infrared spectral range [31]. This remarkable 
property enables the possibility to use the surface signal 
measured at long wavelengths, where atmospheric scatter- 
ing is negligible, to model the surface contribution also at 


(K r K)f t = I< r R s . (9) 

Each reflectance measured by the sensor is instead a 
certain non-linear function of the kernel coefficients. This 
non-linear function is represented by the entire radiative 
transfer model F, and to minimize the mismatch with the 
observations, the residuals <p t 

(fit t) = R i -F i (f) (10) 

are used to calculate the scalar cost function Q as the sum 
of their squares 

<P(f)= ^^.(f) 2 = ^ r (f)^(f). (11) 

i = 1 
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80 

Sf) 


N 


2 Y. < t > i 


8<t>i 

dfj 


= 0 


( 12 ) 


and an iterative approach must be used to solve these 
equations. The problem, written in compact form as 
= argmin f {<P(f)}, may be formidable to solve when 
looking for an absolute minimum of the cost function. 
Several approaches are available when restricting the 
search to a local minimum within a small region <5, so as 
to have O(V) < <P(f)Vf : Ilf— f* II < <5. 

Pending differentiability and good behavior of F and its 
derivatives within the neighborhood of the current itera- 
tion state f\ we can use a Taylor expansion to write 

FA f ) ~ FA f k ) + + °( ll(5fk 1 12 ) < 1 3 ) 

j 

where Sf k =fj-f k and /j- is the Jacobian, or array of partial 
derivatives of F, with respect to each fj. The residuals are 
then approximated by 

M 

~ [Si-FA f fc )J- £ pfj ( 1 4) 

j = i 

from which we also derive dtjf/dfj = -j|. Substituting in 
(12) we obtain the linear set of equations [R— F(f fe )J = J k,r <5f* 
that can be cast in the form of normal equations as in (9): 

(J kJ J k )df k =J fc [R-F(f k )] (15) 

with the difference that now the unknown is the increment 
<5f\ the regressor is the departure from the model R-F and 
the matrix involved in the inversion is the Jacobian. In other 
words, whenever a first-order Taylor expansion is a good 
approximation, the search for a minimum of a non-linear 
function can be turned into a (iterative) linear regression on 
its increments. The advantage of this approach becomes 
clear when considering that the Jacobians are a byproduct of 
the Doubling Adding simulations, conveniently output free 


of additional computational cost. At each iteration step k, a 
better approximation to f* is obtained from the recursive 
relation f k 1 1 = f k + <5f\ and the model run again to obtain a 
new value for Sf k . This procedure is repeated until conver- 
gence, essentially applying Gauss-Newton’s method for the 
search of a local minimum. This method is extremely 
powerful and performs best with smooth functions of low 
curvature. The convergence is quadratic for functions that 
are twice differentiable in the neighborhood of a simple 
root, and can be reached in one iteration in the limit case of 
linear functions. 

The concepts illustrated above are the basis for the 
linearization of radiative transfer codes, obtained for 
example for the Discrete Ordinate [37,38] and the Doubling 
Adding [39] formulations. Specific applications focus on 
trace gas retrievals [39], surface properties [40,41] and 
extensions from plane parallel to spherical atmospheres 
[42,43]. A few recent studies, albeit applied to total 
reflectance only, explore regularization of ill-conditioned 
inversions deriving from a limited set of viewing angles 
[44], a situation that highlights the advantages of multi- 
angle measurements like those provided by the RSP. 

To test its performance, the iterative method was first 
applied to a synthetic scene consisting of a purely Ray- 
leigh atmosphere at 410 nm, where scattering is signifi- 
cant and many of the surface features are washed out. The 
solar zenith angle was set at 30°. The surface was 
prepared according to the model in Eqs. (2) and (5) 
by choosing the following set of parameters: j / iso ,/ vol , 
/geo/fresneil = (0.25, 0.19, 0.05, 0.75}. Additional quantities 
are required to specify the geometric kernel [8], repre- 
senting the ratios of the tree to crown height and of the 
crown vertical to horizontal extent. These parameters 
were chosen as the MODIS default of h/b=2 and f>/r=l, 
respectively, indicating spherical crowns of half the 
height of the tree. The resulting surface reflectance, 
unknown in real measurements and which we ultimately 
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Fig. 1. Gauss-Newton’s iterative search for the parameters optimally estimating the synthetic scene described in the text. The plot on the left is relative 
to the total reflectance, modeled as in Eq. (2); the plot on the right to the polarized reflectance (see Eq. (5)). From knowledge of the reflectance measured 
by the RSP (blue curve), the method converges fast from the initial guess (path radiance-corrected signal, solid cyan curve fitted by the cyan dashed 
curve) to the unknown surface signal (thick solid black curve). Fits during different iterations are represented with dashed red curves, with the converged 
answer of increased thickness. Both for total and polarized reflectance, the residuals get very small after the second iteration already, making it difficult 
to distinguish intermediate fits from the values found in the final iteration. (For interpretation of the references to color in this figure legend, the reader is 
referred to the web version of this article.) 
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want to retrieve, is shown with a black curve in Fig. 1 for 
the total (left panel) and polarized reflectance (right 
panel). The signal simulated for an instrument like the 
RSP, flying at an altitude of 8.5 km along the principal 
plane, is shown by the blue curves. 

An initial guess for the model parameters is needed to 
initialize the iteration. A first estimate for the surface 
(total and polarized) reflectance is obtained by taking the 
difference between the observed signal and the (total and 
polarized) path radiance, defined as the radiance scattered 
into the instrument’s field of view without having inter- 
acted with the surface. A common expedient used in 
radiative transfer to estimate the path radiance is to 
repeat the simulation that generated a scene with a 
surface of zero albedo imposed as a lower boundary 
condition. The total signal corrected for the path radiance 
is shown with solid cyan lines. This approximation to the 
surface signal does not take into account possible effects 
due to multiple interactions between the surface and the 
atmosphere. In this respect, it is worthwhile noting the 
vanishing surface contribution to the polarized reflec- 
tance for directions close to backscatter where Rayleigh 
processes take over. The cusp found for the total reflec- 
tance is instead the clear hotspot signature typical of 
vegetation models. Nonetheless, a first fit with Eqs. (9) 
and (5), shown by the dashed cyan curves, provides 
parameter values close to the real solution enough to 
ensure convergence. This first fit is obtained with a 
Levenberg-Marquardt type of routine, based on the MIN- 
PACK FORTRAN package (http://www.netlib.org/minpack/ 
and translated for the Interactive Data Language (IDL) by 
Markwardt [45], 

When the model is run again with this choice of 
modeling parameters, the surface signal is found to resem- 
ble the first fit (see the almost overlapping red dashed line). 
The Jacobians output by the code are then used in Eq. (15), 
implemented with the IDL LA_LEAST_SQUARES procedure, 
to find the steps <5f that are then used to update the 
parameters. A second iteration produces a significantly 
different signal that already fits the sought surface reflec- 
tance very well. Further iterations account for finer refine- 
ments, barely distinguishable in the plots, until the 
algorithm stops based on a condition imposed on the 
difference in residuals from the previous iteration (not 
greater than 1 0 7 in this specific case, both for total and 
polarized reflectance). Based on the discussion above, this 
fast convergence is an indicator of the low curvature of the 
hypersurface represented by the cost function (10) in 
“residual” space. Given the different order of magnitude 
of the total and polarized reflectance, it is easy to reinforce 


Table 1 

Parameter evolution during optimization of the synthetic scene. 


Step 

/,„ (0.25) 

(0.19) 

/ Ee „ (0.05) 

/, res (0.75) 

GUESS 

0.165154 

0.050942 

0.023309 

0.625167 

ITER 1 

0.248204 

0.193857 

0.048905 

0.950605 

ITER 2 

0.249348 

0.187381 

0.050515 

0.765108 

ITER 3 

0.249909 

0.190033 

0.050013 

0.738954 

ITER 4 

0.250042 

0.190115 

0.049978 

0.750050 

ITER 5 

0.250001 

0.189992 

0.050001 

0.750467 


the convergence criteria by thresholding also the absolute 
values of the residuals. 

Table 1 summarizes the evolution of parameters dur- 
ing the retrieval. 


3. Instrumentation and research flight 


The reflectance terminology involved in remote sen- 
sing can be sometimes confusing and we have collected a 
compendium of relevant quantities in Appendix A. Details 
on the Research Scanning Polarimeter (RSP) are given in 
Appendix B. Here, it suffices to specify that this multi- 
spectral instrument collects multi-angle measurements of 
overflown targets in nine selected bands spanning the 
410-2264 nm wavelength range. For each angle and 
wavelength, the quantities being measured are the first 
three components of the Stokes vector (I, Qand U), which 
are related to reflectance (more precisely to the Bidirec- 
tional Reflectance Factor, BRF) as follows: 


Ri = 


nrf 


F 0 cos i 


(16) 


Ro = 


nrf, 


F o cos e s 


(17) 


Ru = 


nrf 


F 0 cos 0, 


U 


(18) 


where F 0 is the exoatmospheric average solar incident 
flux, r 0 the Sun-Earth distance, and 9 S the solar 
zenith angle. 

During the 2010 Carbonaceous Aerosol Radiative Effect 
Study (CARES) campaign [46] in Sacramento, CA, favor- 
able conditions with snow persisting late in the season 
provided the opportunity to include a few flight legs over 
the Sierra Nevada mountain range in one of the research 
flight of the NASA King Air B-200 on which the RSP was 
deployed. The flight took place on 10 June, and resulted in 
a limited dataset useful to analyze and characterize the 
signal from snowy surfaces (see Fig. 2). Of particular 
interest is the question of whether the polarized reflec- 
tance of snow is spectrally flat as that of other land 
surfaces, because this property would allow to apply the 
same algorithms employed for aerosol retrievals [47], in 
particular optimal estimations based on inverse methods. 

The successful collection of the dataset depended 
both on careful planning and adjustments during flight 
operations. The snow cover had been continuously mon- 
itored via the data downloaded from the NOAA 
National Operational Hydrologic Remote Sensing Center 
(http://www.nohrsc.nws.gov/interactive/html/map.html). 
Although it was not possible to collect snow samples, a 
few days prior the flight areas with similar snow condi- 
tions in the South Tahoe region (Ralston Peak, Eldorado 
National Forest, CA) were surveilled on foot and the snow 
was found granular and with considerable water content, 
as shown in Fig. 2. 

The flight legs were designed pondering the best 
tradeoff among sufficient snow cover, flatness of terrain 
and orientation of the principal plane: flying in this 
direction maximizes the range of scattering angles 
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Fig. 2. Flight trajectory (shaded vertical curtain) leading to the collection of the polarimetric dataset over snow. North is at bottom left. Coming from 
Mono Lake (upper left), we overflew the high country past Mt. Warren onward to North Peak on the Sierra Divide (light green), which at this location also 
marks the border between Inyo National Forest and Yosemite National Park. The real-time panorama swath was created by blanketing the 3D terrain in 
Google Earth with a merge of images captured by the B-200 downward looking camera. The inset shows the appearance of snow on Ralston Peak 
(Eldorado National Forest, CA) on 6 June, with wet clusters of nearly spherical grains a few millimeters in diameter. Although the sampling location was 
only relatively close in space and time, these late-season snow conditions were likely to be similar to those scanned during the research flight. (For 
interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 


spanned during RSP scans. Also, wind forecast were 
monitored since head or tail winds limit aircraft yaw, 
which prevents features limited in size from being 
scanned at all available angles. Experience acquired dur- 
ing the analysis of RSP data from previous campaigns 
confirms that corrections for “crabbing” are indeed a 
major challenge during data aggregation. 

Airborne determination of the reflection properties of a 
surface benefits from low-altitude passes, so as to max- 
imize resolution and minimize atmospheric interference. 
Nevertheless, stable aircraft course is challenged by the 
turbulence associated with a rugged topography. With 
time constraints preventing the possibility of repeated 
overpasses, we chose to fly at an altitude that preserved 
the High Spectral Resolution Lidar (HSRL [48]) operations, 
since this instrument provides vertical characterization of 
the atmosphere below the aircraft but cannot be operated 
below altitudes established by Federal Aviation Adminis- 
tration regulations on eye-safety. 

The most interesting data were collected along a flight 
segment from Mono Lake to Yosemite National Park (see 
Fig. 2 and Appendix C). Curiously, this area belongs to 
regions surveyed in the limited set of published works 
that were found to deal with hyperspectral imaging of 
alpine snow in the United States [49]. Regional variability 
and interconnection between grain size and wetness were 
investigated based on Airborne Visible InfraRed Imaging 
Spectrometer (AV1RIS) data collected near Tioga Pass [49], 
and the effect on spectral reflectance of red algae 
( Chlamydomonas nivalis ) blooming onto old snowfields 
was observed by the same instrument on the drainages 
of Mt. Conness [50,51], 

At the time the transect took place, the air was 
exceptionally clean of aerosols with an aerosol optical 


depth below the aircraft measured by the HSRL equal to 

0. 015. (at 532 nm). Last-minute, in-flight adjustments 

worth mentioning included: 

1. descending at lower altitude, in this case the minimum 
allowed for HSRL operations ( ~ 7.5 km). Flying over 
plateaus with an average elevation of 3200 m, the RSP 
resolution was ~ 60 m. From the same altitude, scanning 
a target throughout the whole RSP field of view ( + 60 ) 
required a straight leg of a length x = 2h tan 60 ~ 7 km 
(5.3 km for the reasonable range of +40°); 

2. disengaging the autopilot, which tends to overshoot 
the corrections to variations in attitude. This choice 
was particularly fortunate as manual flight greatly 
stabilized the orientation of the scanning plane 
throughout subsequent scans; 

3. reducing the speed, which also improves resolution. 
This maneuver was beneficial only due to the flight leg 
being aligned with the wind direction; 

4. deviating approximately 20° from the planned course, 
to overfly the flattest terrain. Varying the course 
according to visually identified features is difficult 
since the visibility from the cockpit is limited in the 
downward direction. Also, the variation must be esti- 
mated several kilometers in advance to make sure the 
feature of interest is overflown within a straight leg. 


4. Optical properties of snow 

Snow can be physically considered as a lattice of ice 
grains, whose radii are typically between 50 pm and 
1000 pm depending on age, although significant amount 
of liquid water can be present especially after the onset of 
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melting. The complex indices of refraction of water and 
ice, together with the size and shape distribution of the 
ice grains, therefore determine the optical properties of 
this highly scattering medium. The overall picture is 
complicated by the assortment of grain types [52,53] 
and a remarkable tendency to metamorphism, at the 
point that characterizing or even defining size and shape 
is often a struggle for both measurements and models. 
Moreover, topographic effects and atmospheric agents 
can easily alter the essentially isotropic properties of 
fresh snow resting on a flat ground. Snowfields deposited 
on rugged terrain, typically associated with alpine envir- 
onments, receive special attention by the remote sensing 
community motivated by the need of inferring important 
information on the hydrological cycle [54,49]. Such sce- 
nario can be challenging to model, because radiative 
transfer codes typically rely on the assumption of a 
bottom boundary described by a surface located at a well 
defined altitude and of unlimited horizontal extent. 
Further caveats are represented by sub-pixels heteroge- 
neity [49,55] and shadowing effects due exposed soil and 
vegetation. 

For wavelengths in the visible region of the electro- 
magnetic spectrum, snow is the natural surface of highest 
reflectance, with a spectral albedo that can reach values 
up to 0.98 for fresh dry snow [56], The white appearance 
is due to nearly constant behavior of the real part 
of the refractive index, which is also very similar for both 
the ice and water phase with that of the ice being only 
slightly shifted towards longer wavelengths. The absorp- 
tion coefficient (imaginary part of the refractive index) is 
also similar, and more importantly is very small prevent- 
ing liquid water from appreciably affecting the 
reflectance. 

As the spectral window shifts towards the near-infra- 
red and up to 2.5 pm, the properties drastically change 
and the “color” of snow becomes all but monotone. The 
absorption coefficients of both ice and water experience a 
similar, remarkable increase covering seven orders of 
magnitude. As a consequence, the reflectance plummets, 
to partially recover only around 1.9 and 2.4 pm. Very dark 
intervals are found at 1.5 pm and 2.0 pm. Differently than 
in the visible, the presence of liquid water can affect the 
reflectance, also outside the 1.35-1.75 pm region where 
ice is slightly more absorptive. The reason is that liquid 
water promotes clustering, leading to larger ice grains 
(up to a fraction of a centimeter) of complex shape. 

Snow albedo increases with the solar zenith angle. This 
effect is detectable in the visible, and very evident in the 
near-infrared, where the strongly forward-scattering 
ice grains (especially those of larger size) have less of 
a chance to attenuate via absorption a beam incident 
at a grazing geometry than they have to attenuate a more 
penetrating normal beam. The higher probability of 
absorption determined by longer in-grain paths also 
explains the sensitivity to grain size, which is greatest at 
1.0— 1.3 pm and has been successfully modeled and 
applied to observational datasets [57-59]. All algorithms 
exploit the ratios between relatively unaffected bands in 
the visible and largely affected bands in the Near Infra- 
Red (NIR). 


In terms of angular behavior, the most prominent 
observed feature is the increase in reflectance at forward 
scattering angles. There is some indication that this 
increase is partially suppressed by both surface roughness 
and wetness, but accurate quantification and accepted 
explanation for these effects is still lacking [60], Different 
water phases can be detected exploiting the signal col- 
lected in water absorption bands by hyperspectral instru- 
ments [61]. This task is of particular interest, because of 
the importance of monitoring snow melting processes. 
A weaker increase in reflectance is measurable also in the 
backscatter region, and can be enhanced by roughness. 
Specular signatures are missing from most of snow 
samples analyzed in the published literature [60], In some 
cases, wet samples show weak structures [62], but peaks 
in the specular direction were only detected in rough ice 
layers, for which the mirror effect of macro-facets is the 
most logical explanation. 

These arguments obviously hold for a pure (at most 
contaminated by water) and well formed (semi-infinite) 
snowpack. Snow depth and impurities deposited over or 
embedded within the snow layer have a large effect on 
the reflectance: the former because of partial exposure 
of the underlying surface, the latter depending on the 
radically different refractive index of black particulate 
such as soot. The thickness required to make a snow layer 
optically thick depends on the grain size and the wave- 
length. A layer of a few millimeters is sufficient to mask 
the underlying surface at long wavelengths, while a few 
centimeters of fresh snow to tens of centimeters of wet 
snow are required to fulfill the same approximation in the 
visible. The effect of impurities can at times be drastic. 
While the high absorption properties of ice in the NIR are 
little affected by the presence of extraneous absorbers, in 
the visible even a few parts per billion by weight of dark 
particles are sufficient to contaminate the transparency 
and lower the albedo. This phenomenon has received 
much attention due to its involvement within a feedback 
mechanism leading to increased melt rate which can 
affect climate on a global scale [63,64], 

Radiative transfer simulations exploring the sensitivity 
to different parameters were prompted by the seminal 
studies of Warren and Wiscombe [65,66], and continued 
in the following decades [67-73,57], Very few studies 
extended the modeling to the full Stokes vector [74,75], 
also because the scarcity of observational datasets 
[76-78], Most of the knowledge cumulated so far comes 
from in situ measurements that exploit goniometers in 
the field or in the laboratory [78,79]. The degree of linear 
polarization was in any case measurable with satisfactory 
signal to noise ratio only in the forward scattering region, 
and shows to be weak in the visible, getting stronger in 
the infrared. This fact is explained by the fewer scattering 
events undergone by long wavelengths before they are 
reflected off the surface; at visible wavelengths larger 
penetration depths allow multiple scattering to suppress 
polarization. Specular reflection peaks are associated to 
strong polarization, but as previously mentioned they 
have only been detected in frozen, macroscopically rough 
layers (and only partially in wet samples). More investi- 
gations are needed to establish the weight of this claim, 
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especially regarding the degree of polarization measured 
around the Brewster’s geometry (roughly 52° of incidence 
for ice), for which polarization by reflection is expected to 
be maximized. On the other hand, a most interesting 
feature that seems confirmed is the significant reduction 
of the forward polarization by wet snow (as opposed to 
dry old snow), which could prove diagnostic for the 
detection of water content. 


5. Data analysis and discussion 


The negligible aerosol optical depth measured by the 
HSRL (see Section 3) greatly simplified our analysis since 
the atmospheric correction consists only of the removal of 
the Rayleigh contribution. In any case, accurate analyses 
of RSP data rely on correction for band-specific gaseous 
absorption. The radiance measured in the Short-Wave 
InfraRed (SWIR) channels at 1594 and 2264 nm need be 
corrected for water vapor absorption. The ratio between 
the signal in a channel subject to water vapor absorption 
and one in a window can be used over bright areas 
to retrieve the precipitable amount [80,81]. The RSP offers 
this possibility through the reflectances measured at 
960 nm and 864 nm which were used to calculate 


R,(960) _ i(amw/ 
R4864) P 


(19) 


where AM is the airmass, W is the water vapor column (in 
precipitable cm), and a and /? are instrument-specific 
constants [82], 

The possibility of utilizing to this purpose the bright 
scans over snow was enticing. Nevertheless, an ice/water 
mix typical of late-season snow can spoil the usefulness of 
the ratio, because the imaginary part of the refractive 
index of ice is appreciably less than that of water at 
960 nm (3.5xl0“ 6 versus 7.5xl0“ 6 [83,84]). The 
amount of precipitable water retrieved with the proce- 
dure was 1.16 cm, indeed an implausible value for the 
4 km of atmosphere between the aircraft and the ground. 
Repeating the procedure for scans over the Tioga Crest 
(mostly exposed, talus/scree slopes of metamorphic rocks 
at 3600 m asl, see Appendix C) and over sunglint reflected 
off Mono Lake (1950 m asl), the retrieved value was 
0.15 cm. This low value allowed to skip the correction 
for water vapor for the two RSP bands. 

Minimal corrections for N0 2 and 0 3 , amounting to less 
than 0.5%, were applied the 410 nm and 469 nm channels, 
based on the concentration of these gases retrieved form 
the Ozone Monitoring Instrument (OMI [85]) and TOMS. 

It is recommended [86,9] that if the inversion produces 
a negative weight for one of the total reflectance kernels, 
the inversion be repeated with that value constrained to 
zero. This sequence ensures that the quasi-orthogonality 
of the kernels is fully exploited, minimizing the risk of 
unstable inversion and/or retrieval of unphysical reflec- 
tances. Since the analytical form of the Jacobians output 
by the DA code corresponds to dR(f)/d log f k , the solution 
to Eq. (15) actually provides Slogf. To constrain the 
search to positive values of the kernels’ weights, this 
step size can be used to update the parameters within a 


logarithmic-exponential transformation 
exp[log f k + ^ log f k ] = exp[log f k + <5f k /f k ] 

= f fc exp[<5f k /f k ] ~ f ,( [l + <5f k /f k ] = f k + 5f k (20) 

A discussion on the applicability of this kind of 
transformation to parameter estimation is given in 
Dubovik and King [87], 

Some authors [A. Lyapustin, personal communication ] 
prefer instead to regard the kernels as purely mathema- 
tical functions and accept negative weights; the ideal 
treatment is still a source of debate. The geometric and 
volumetric kernels were not designed for a peculiar sur- 
face as is the case for snow, whose reflectance is largely 
dominated by the isotropic component. Small deviations, 
occasionally observed at very large angles and mostly 
affected by scanning uncertainties on a rough terrain, can 
introduce large variability in the weights, including small 
negative values indeed observed for some of the channels 
in our analysis. A complete characterization of snow total 
reflectance model is hampered by the unique nature of 
this limited dataset. Since the primary goal of this study is 
to quantify the spectral variability of the polarized reflec- 
tance, negative occurrences were accepted as were, and a 
thorough study of the performance of surface reflectance 
models on snow was delayed to a more representative 
dataset. 

Once the iterative procedure has converged, the diffuse 
component calculated by the model can be subtracted 
from the observations, and the difference divided by the 
two-pass transmittance to obtain a best estimate for the 
surface reflectance. This result overlaps to the converged 
answer if the model is capable of an accurate description 
of the reflection of the diffuse component throughout 
the range of viewing angles. As observed for other land 
surfaces [47,88,89], the Fresnel model performed worst 
close to the backscatter direction, where both the atmo- 
sphere and the surface introduce very small polarization 
often characterized by inversions. In this angular region, 
due to the virtually identical reflectivity predicted at 
backscatter (perpendicular incidence) for the parallel 
and perpendicular components, the theoretical polarized 
reflectance instead decreases monotonously to zero. 

In the following we contrast the spectral behavior of 
the snow surface polarized reflectance obtained after 
convergence of the iterative scheme at each wavelength 
with that of a standard atmospheric correction. 

The standard scheme is based on subtraction of direct 
and diffuse atmospheric contributions, including an esti- 
mate of the multiple interactions coupling the surface and 
the atmosphere. For the polarized reflectance, we tested 
an ad hoc approach used in past surface polarized reflec- 
tance retrievals from POLDER and RSP measurements 
[13,90]. Neglecting multiple interactions, the polarized 
reflectance R P (A,z) measured by a polarimeter at an 
altitude z can be written as 

R P U,z) = R“ tm (2,z)+T u (x,z)Rp r/ a) (21) 

where R P tm (2) and R/^U) are the path and surface polar- 
ized reflectance, respectively (see Section 2). The latter is 
weighted by the two-pass transmittance from the top of 
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the atmosphere to the surface and up to instrument 
altitude for each viewing angle, T u (x,z). Note that we 
have maintained the vectorial notation to remember that 
we are dealing with arrays of measurements at different 
viewing geometries (scattering angles), and reintroduced 
the explicit dependence on wavelength. In absence of 
appreciable aerosols, the transmittance is determined by 
the molecular optical depth t m and the coupling effects 
can be approximated by rescaling T m with a factor i/c 
From Eq. (21), i/r is obtained through a linear regression on 


R P (l,z)-Rp tm (l,z) = Rp^(l)exp 




( 22 ) 


viz., 


<A = 


l 

logT n (l,z) 


log 


RpU,z)-R“ tm (/l,z) 


(23) 


Attempts to reproduce the factors reported in the 
literature [13,90] found that the latter were consistently 
overestimated. A possible reason for this lack of reprodu- 
cibility is that those studies retrieved the factor t/r in 
presence of significant aerosol loads, which dominate the 
polarization signatures therefore impacting the accuracy 
with which the corrective term relative to the Rayleigh 
contribution is determined. 

The retrieved total surface reflectance decisively points 
to the typical signatures of snow, with high intensity 
measured in all three scenes by the visible channels and 


the expected drastic decrease in the SWIR. More impor- 
tantly, the variations in visible “albedo" are remarkable 
among the three scenes. It has to be noted that what 
plotted in Fig. 4 are rigorously speaking apparent reflec- 
tances. In particular, the western face of North Peak 
(Scene C) is characterized by a constant slope of 
about 40°. Given the close alignment of the direction of 
steepest slope to both flight trajectory and principal 
plane, a simple correction to the illumination geometry 
can be introduced by adjusting the solar zenith angle in 
Eq. (B.l). This correction for Scene C can amount to a 80% 
reduction of the reflectances and bring it closer to Scenes 
A and B. Scene B presents a slightly upward slope as well, 
which once corrected for would increase the reflectance. 
It is nevertheless unlikely that the apparent character of 
these reflectances is the sole responsible factor for the 
observed variability, especially considering the differ- 
ences assessable by the naked eye among to the three 
snowfields. 

It is noteworthy that Scene B (Fig. 3) exhibits the 
highest reflectance at 670 nm, making the visible color of 
this snowfield slightly pink. This behavior can possibly be 
attributed to the presence of the snow algae common in 
this area, as discussed at the end of Section 3. 

Finally, the presence of water/ice mix agrees with the 
enhanced backscatter measured over the analyzed snow- 
fields, although a major contribution to this feature can 
be also ascribed to surface roughness as observed by 
Peltoniemi et al. [60]. 


Scene A (Warren Mt.) 


Scene B (20 Lakes Basin) 


Scene C (North Peak) 


Scattering angle (deg) 

172 164 154 144 134 124 115 105 


Scattering angle (deg) 

172 164 154 144 134 124 114 104 


Scattering angle (deg) 

172 163 154 144 134 124 114 104 




-30 -20 -10 0 10 20 30 40 

Viewing angle (deg) 



-30 -20 -10 0 10 20 30 40 

Viewing angle (deg) 


-30 -20 -10 0 10 20 30 40 

Viewing angle (deg) 


LU 

o 


o 

LU 


LU 

IT 


£ 

o 



Viewing angle (deg) 


Viewing angle (deg) 


Viewing angle (deg) 


Fig. 3. Surface reflectance retrieved with the iterative procedure described in Section 2. The three columns are relative to the three selected scenes 
(cf. Fig. 2). The top row pertains to the polarized reflectance, and the bottom row to the polarized reflectance. Note the large variations in total reflectance 
among the three scenes. The spectral variability and the absolute value of the polarized reflectance are instead small. 
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Fig. 4. The residuals in the top row of this figure correspond to the differences between the surface polarized reflectances retrieved for each RSP 
wavelength (upper row of Fig. 3) relative to the 2264 nm channel. Away from the backscattering region, where land models are known to fail, similar 
results are obtained via a standard correction based on the simplistic subtraction of atmospheric scattering contributions (bottom row). 


For the polarized reflectance, the two methods give 
very similar results, with residuals limited to 2xl0~ 3 
except for the 670 nm and 864 nm in Scenes A and B. This 
value is 5-10 times higher than that found for other land 
surfaces. It can also be seen that the model (upper row) 
does not adequately describe the behavior in the back- 
scatter region, i.e. at scattering angles > 160°. 

Recognizing the similarity between the angular beha- 
vior of the polarized reflectance retrieved over the sur- 
veyed snow and that of rough ice clouds, we extended our 
study attempting to retrieve surface descriptive para- 
meters. The approach is analogous to that described in 
a recent publication where RSP measurements were 
exploited for the retrieval of optically relevant parameters 
of ice crystals forming cirrus clouds [81], Firstly, a 
database of optical properties of single plates and col- 
umns was constructed based on geometric optics [91] 
with size, aspect ratios and small-scale roughness as free 
parameters [92-94]. Then radiative transfer simulations 
were run for a layer of ice particles of high optical depth 
(200) occupying the first meter of the model atmosphere. 
The departure from the data of the polarized reflectance 
simulated for each habit in the database (Root Mean 
Square Error, RMSE) was stored, and used as a measure 
of the goodness of the fit. A sample result of this 
procedure is presented in Fig. 5 for Scene B. Here we 


chose to exploit the polarized reflectance at 864 nm to 
minimize uncertainties associated with absorption esti- 
mates, since ice does not appreciably absorb at this 
wavelength. The upper panel of Fig. 5 shows that good 
fits are obtained for a variety of ice models. The lower 
panel shows a contour plot of the RMSE. Looking for 
minima (marked with red diamonds), it is observed that 
all optimal fits pertain to crystals affected by a surface 
roughness of 0.3-0.35. The aspect ratios are instead not 
unambiguously retrieved and include a collection of 
extreme geometries (i.e., very long columns), a fact that 
can potentially be exploited when attempting the con- 
struction of a surface model. Large departures occur from 
the model curves relative to pristine or nearly-pristine ice 
crystals, since these habits generate well-known halos not 
present in the analyzed data. 

A final argument is especially important to put this 
work in perspective, and it is summarized in Fig. 6. Here 
we show with solid lines the spectral behavior of the 
polarized reflectance acquired by the RSP just 20 min 
after the B200 surveyed the snowfields. At this time, the 
aircraft was overflying the Fresno area where the HSRL 
measured an aerosol optical depth close to 0.1 at 532 nm. 
The dashed lines result from multiplying the surface 
signals retrieved for the snow (Scene B) by the transmit- 
tance factor expf-T^/^o-F/O], and therefore approximate 
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Fig. 5. Search for a best fit to the polarized reflectance measured over Scene B. The snow layer was treated in radiative transfer simulations as a thick ice 
cloud at the surface, and the ice crystals descriptive parameters let vary according to a database calculated via geometric optics. The lower panel shows a 
contour plot of the RMSE (departure from the model) as a function of model parameters. All best fits (red diamonds) correspond to grains characterized 
by a roughness of ~0.3 and not well-determined, extreme aspect ratios. (For interpretation of the references to color in this figure legend, the reader is 
referred to the web version of this article.) 


the surface contribution as observed from the RSP altitude 
if the ground around Fresno had been covered in snow. 
The spectral AOD was determined for a simple aerosol 
model (r eff = 0.15 pm, v eff = 0.15), with a particle number 
concentration adjusted so as to yield the HSRL-measured 
optical depth. It should be noted that as long as the single 
scattering approximation holds, the polarized reflectance 
is proportional to the optical depth itself (see, for exam- 
ple, Eq. (2.1) in [16]). For this moderate aerosol load 
already, and more so at increasing optical depths, the 
surface contribution to the measured polarized reflec- 
tance is small over a large portion of the angular range. 
With the exception of the 670 nm and 864 nm bands, at 


which the analyzed snowfield exhibited somewhat pecu- 
liar behavior, the spectral variability is also limited. This 
fact suggests that, pending a more thorough characteriza- 
tion based on data collected over flatter snow surfaces, 
polarization measurements can be exploited to improve 
the estimate of the surface contribution, thereby increas- 
ing the accuracy of aerosol retrievals. 

6. Conclusions 

Airborne polarimetric data over snow, with the accu- 
racy provided by an instrument such as the RSP, are to our 
knowledge unprecedented. This paper reports the first 
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Fig. 6. Comparison between the small magnitude of the snow polarized 
reflectance (dashed lines) and the polarized reflectance measured 
20 min later in the flight (solid lines) over the Fresno area, characterized 
by an AOD 532 = 0.1. The RSP 1590 nm band is excluded for clarity. 

known investigation of the behavior of the polarized 
reflectance of snow over the 410-2264 nm spectral range, 
of interest to both the cryosphere and aerosol remote 
sensing communities. The dataset was collected during an 
RSP survey over alpine snowfields in the Sierra Nevada 
range, in June 2010. The surface signal retrieved via a 
standard atmospheric correction is compared with the 
results obtained from an iterative correction scheme, 
which automatically isolates the surface contribution 
from remote sensing polarimetric observations. Based on 
the linearization of a Doubling-Adding radiative transfer 
code with respect to the surface descriptive parameters, 
this method eliminates the need of an atmospheric 
correction. 

The primary conclusion of this study is that the 
polarized reflectance of the analyzed snow type is not as 
spectrally neutral (“gray”) as that of other land surfaces, 
but the differences are nonetheless small. The absolute 
value of polarized reflectance is also small. These two 
facts together augment the confidence in the usefulness of 
polarimetric measurements to effectively separate surface 
and atmospheric contribution in the analysis of remote 
sensing measurements. The application of the iterative 
correction scheme highlights how current surface models 
do not fully reproduce the polarization signatures of 
snow, especially near backscatter where the models are 
known to fail even for the traditional land surfaces for 
which they were designed. 

We have also attempted to model the snow layer in 
radiative transfer computations as an ensemble of ice 
crystals at the surface. The results show that the polariza- 
tion signatures of the analyzed snow fields resemble 
those of an optically thick, rough ice cloud. 

The acquisition of the dataset was a particularly 
valuable experience because of the challenges presented 
when scanning spatially limited targets on uneven ter- 
rain, as opposed to horizontally homogenous aerosol 
layers or cloud decks. Successful data collection was 
possible due to over 10 years of experience in RSP 
research flights, and is auspicious for the collection of a 
set of dedicated measurements. At present, the RSP is 
arguably the best available asset for high-accuracy, air- 
borne characterization of the polarized reflectance of 


snow and ice landscapes, and consequently for the deter- 
mination of optical and microphysical descriptive para- 
meters of both the surface and the aerosols above it. 
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Appendix A. Radiometric definitions 

This appendix lists a series of definitions pertinent to 
reflectance measurements consistent with the nomencla- 
ture recommended by Nicodemus et al. [95], albeit 
tailored to remote-sensing applications as in Schaep- 
man-Strub et al. [96] and Thomas and Stamnes [97]. 

Electromagnetic radiation transports energy through 
space at a rate determined by the flux of the associated 
Poynting vector [29]. At any point in space and time, light 
can therefore be regarded as radiant energy, Q, measured 
in Joules (J). A spectral counterpart can be defined for 
energy and all the quantities discussed henceforth, by 
considering the derivative with respect to wavelength. 
This dependence is dropped for simplicity from the 
notation, since it is straightforwardly reintroduced when- 
ever the application requires it. 

The energy per unit time is the radiant flux or power, 
and has the units of Watts 

$= d £m. (A.1) 

The Radiant flux density is the radiant flux per unit area 

E=^(Wm- 2 ) (A.2) 

and is also called irradiance, or radiant exitance M when 
considering radiation originating from a surface. The 
radiance L embodies the idea of a pencil-beam of radiation 
and is defined as the infinitesimal amount of radiant flux 
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contained within a solid angle da> oriented in a specified 
direction (0 r , <j> r ), and crossing a surface dA: 


L = 


d 2 0 


cos 9 dA doj 


( W m 2 sr 1 ). 


(A.3) 


Here cos 8 is the angle between the surface normal and 
the direction of propagation, so that dA cos 9 represents 
the area projected onto a plane normal to the “beam". 
This definition implies that the radiance is independent 
on the particular orientation of the surface dA. It also 
implies that the irradiance is obtained by integrating the 
radiance over all the hemisphere: 


E = 



L(6 r ,(j) r )sin 8 r cos 8 r d8 r d(j> r 


(A.4) 


which in the specific case of isotropic illumination 
reduces to E = nL. 

Radiance is sometimes misleadingly referred to as 
intensity. The ( radiant ) intensity I is instead correctly 
defined as the amount of radiant flux contained in the 
elemental cone: 

d<P , 

/= — (Wsr- 1 ). (A.5) 


To describe the local reflection properties of a surface, 
the light scattered by a surface element in the direction 
(0 r , 4> r ) must be related to that incident on the same 
surface element from the direction (0>, </>,). The radiance 
traveling in an infinitesimal cone is dl = L dco. The ratio of 
the radiance reflected by a surface element to the power 
(irradiance) which invests it is the Bidirectional Reflection 
Distribution Function (BRDF): 


BRDF = 


dL r 

L t cos 9 t dcQj 


dL r 

dEi 


(sr^ 1 ). 


(A.6) 


This theoretical definition contains undetectable infi- 
nitesimal quantities, and is therefore not applicable to real 
measurements. Photometers deployed for atmospheric 
studies measure a radiant flux that has generally under- 
gone scattering events: the true physical quantity being 
measured is thus the fraction of incident flux (irradiance) 
reflected in the direction of observation (exitance), or the 
reflectance. A common way to express a reflectance mea- 
surement is through comparison with the response that an 
ideal Lambertian surface with a spherical albedo equal to 1 
would give under the same illumination condition. The 
ratio between the two is the Bidirectional Reflectance 
Factor (BRF), which coincides with the reflectance mea- 
sured by the RSP (see Eqs. (B.1)-(B.3)) since the variations 
in surface properties across the narrow IFOV (for nominal 
flight altitudes) can be assumed to be negligible. 

Several derived quantities are tailored to specific appli- 
cations. The Directional-Hemispherical Reflectance (DHR), 
also known as planetary albedo, is obtained by integrating 
the BRDF of a plane-parallel beam over all viewing angles. 
The Bidirectional-Hemispherical Reflectance (BHR) is also 
called spherical albedo and results from the integration of 
the BRDF over both illumination and viewing angles. In the 
MODIS community, these quantities are referred to as 
“black-sky” or “white-sky” albedo, respectively (the latter 
under isotropic illumination conditions). 


Appendix B. The research scanning polarimeter (RSP) 

The Research Scanning Polarimeter (RSP [14]) was 
built at SpecTIR, under commission of the NASA Goddard 
Institute for Space Studies (GISS) science team, to collect 
data on aerosol and clouds useful for climate studies. It 
has now over 13 years of history of research flights, and 
successful operations prompted its selection as a proto- 
type for the Aerosol Polarimetry Sensor, that was unfor- 
tunately lost as a consequence of a rocket failure during 
the launch of the NASA Glory satellite mission on 4 
March, 2011. 

The RSP optical assembly consists of three pairs of 
boresighted refractive relay telescopes, each pair collect- 
ing measurements in three spectral bands. Of the result- 
ing nine channels, six have centerbands in the visible and 
near infrared portion of the spectrum (410, 469, 555, 670, 
864 and 960 nm) and exploit blue enhanced silicon 
photodiodes, whereas three are in the shortwave infrared 
(1594, 1880 and 2264 nm) and exploit HgCdTe detectors 
cooled to 150 K. 

The paired design is required to determine the I, Qand 
U Stokes parameters from the simultaneous measurement 
of the two orthogonal polarization states at 0° and 90° to 
the meridional plane of the instrument (first telescope), 
and the equivalent intensities at 45° and 135° (second 
telescope). The orthogonal polarization states of the 
incoming light are spatially separated by a Wollaston 
prism and measured simultaneously. This approach 
ensures that the polarization signal is not contaminated 
by uncorrelated spatial or temporal scene intensity varia- 
tions during the course of the polarization measurements, 
which could create “false” polarization. 

The instantaneous field of view (14 mrad) is continu- 
ously scanned within + 60° from nadir by a polarization- 
insensitive scan mirror system, consisting of two mirrors 
each used at 45° angle of incidence and with their planes 
of incidence oriented orthogonally. This design ensures 
that the polarization orientation that is perpendicular to 
the plane of reflection at the first mirror is parallel to 
the plane of reflection at the second mirror so that all 
polarization states are transmitted equally. Scanning 
along-track, the same scene is therefore viewed from 
multiple angles in successive scans. 

In remote sensing applications, it is practical to work 
with reflectance quantities. The first three components of 
the Stokes vector I, Q. and U, measured by the RSP, are 
related to dimensionless Bidirectional Reflectance Factors 
(BRFs, see Appendix A) as follows: 


F 0 cos 9 S 

(B.l) 

< Q 
F 0 cos 9 S ^ 

(B.2) 

nr o u 

F o cos 9 S 

(B.3) 


where 9 S is the solar zenith angle. The annual average 
solar flux incident at the top of the atmosphere, F 0 , is 
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Table B.l 

Absolute error components applying to RSP quantities, with /i 0 indicat- 
ing the cosine of the solar incident angle. The Degree of Linear 
Polarization (DoLP) is defined as R p /R,. 


Quantity 

^inst 

0"cal 

ffpol 

R, 

\fw~ 7 R,/p 0 

0.03 • R, 

- 

R« 

\l 10- 7 R,//r„ 

0.03 • R q 

0.001+0.001 ■ DoLP 

Rp 

\l 10- 7 R,//t 0 

0.03 • R p 

0.001+0.001 ■ DoLP 

DoLP 

/ 10“ 7 

- 

0.001+0.001 ■ DoLP 


V R+„(l+DoLP 2 ) 




corrected for the Sun-Earth distance r 0 (in AU) at the time 
of measurement in order to obtain the proper irradiance. 

Three sources of uncertainty affect the measured 
signal. They are: 

• the instrumental noise, limited by shot-noise and 
depending on R, and incident angle. It is equal for R,. 
R q and the polarized reflectance R p = ^R Q 2 +Ru2 . The 

expression for the Degree of Linear Polarization (DoLP) 
can be found by error propagation; 

• the 3% uncertainty due to absolute radiometric 
calibration-, 

• the uncertainty in the relative polarimetric accuracy, which 
increases with the degree of linear polarization, yet is 
maintained within 0.2% by the RSP internal calibration 
system measuring orthogonal polarization states. 

The choice of polarized reflectance, Q. component, or 
DoLP as a descriptor for the polarization signatures of a 
specific scene depends on the particular scopes of each 
analysis, but the DoLP has the clear advantage of being 
the most accurate measurement since it is not affected by 
the radiometric calibration error. Table B.l reports how 
the uncertainties described above impact each of these 
quantities. 

Appendix C. The Tuolumne Intrusive Suite 

The flight transect analyzed in this paper intersected 
the Tuolumne Intrusive Suite, a geologically highly com- 
plex area at the eastern margin of the Sierra Nevada 
Batholith, intruded into marine sedimentary and volcanic 
rocks [98-100]. The intrusion of the granitic batholith is 
responsible for the development of intense contact meta- 
morphic processes in the country rock. The geologic 
framework reflects in the main topographic features, with 
ridges preferentially aligned NNW-SSE, following the 
major lithologic contacts. Moving E-W from Mono Lake 
to North Peak, Paleozoic ( ~ 500-300 Ma) marine meta- 
morphosed rocks are founded in large, isolated bodies 
(e.g. Mt. Warren), surrounded by granitic and granodiori- 
tic rocks (Gilcrest Peak-Mono Dome alignment). To the W, 
a narrow NNW trending corridor connecting the Tioga 
Crest to Tioga Peak is underlined by late Paleozoic-early 
Mesozoic marine sandstone and limestone, and their 
metamorphic equivalent (marbles). At Saddlebag lake 


the metamorphosed marine rocks give way to a complex 
assemblage of Mesozoic volcanics and volcanoclastics, 
locally highly metamorphosed, that extends up to the 
eastern slope of Shepherd Crest East. From here to the W, 
the granitic-granodioritic lithology of Sierra Nevada Bath- 
olith dominates almost uninterruptedly. The high eleva- 
tion of the area together with a long history of glaciations 
is responsible for the presence of abundant rock debris 
and talus deposits. Bedrock is primarily exposed on steep 
cliffs, crests, and in areas of granitic outcrops. The rest of 
the landscape is covered by scree that can reach several 
meters of thickness in areas of gentle slope or at the base 
of hillsides. Valley bottom of major drainages are draped 
by a variable thickness of glacial deposits (till). 
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